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10.  ABSTRACT  fCk'llna  an  ravaaa#  aMa  II  nacaaaaar  aN  Mwwlfp  Af  Alack  numbme) 

In  recent  years,  the  Kalman  filter  has  been  utilised  extensively  for  passive  target 
motion  analysis  (TMA)  —  an  application  in  which  filter  divergence  is  a  common 
problem.  Available  methods  for  eliminating  divergence  ultimately  involve  increas¬ 
ing  filter  sensitivity  by  discounting  the  influence  of  past  data.  However,  this  pro¬ 
cedure  matos  the  filter  more  susceptible  to  random  errors;  therefore,  to  avoid  un¬ 
necessary  sacrificing  of  noise  performance,  adaptive  control  is  required.  In  this 
report,  the  Kalman  filter  equations  are  derived  and  the  associated  data  measure-  ^ 
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20.  ABSTRACT  (Cont'd) 

control.  An  Important  relationship  between  the  system  performance  index  and 
the  data  residuals  is  established.  By  exploiting  this  relationship,  pertinent 
statistical  properties  of  the  performance  index  are  deduced  and  utilized  as  a 
basis  for  formulating  practical  adaptive  control  criteria.  A  simulation  example 
is  presented  to  demonstrate  divergence  (e.  g. ,  tracking  of  a  maneuvering  tar¬ 
get)  and  significant  improvement  in  performance  is  noted  when  adaptive  control 
is  appended. 
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THE  UTILIZATION  OF  DATA  MEASUREMENT  RESIDUALS 
FOR  ADAPTIVE  KALMAN  FILTERING 

INTRODUCTION 

In  recent  years,  the  Kalman  filter  and  other  recursive  estimation  schemes 
have  been  utilized  extensively  for  target  motion  analysis  (TMA).  Unfortunately, 
in  these  applications,  a  common  problem  known  as  filter  divergence  is  often 
encountered.  Divergence  occurs  when  the  calculated  error  covariance  becomes 
inconsistent  with  the  actual  error  covariance.  Although  there  are  many  possi¬ 
ble  causes  for  divergence,  a  common  source  is  system  modelling  errors.  This 
is  particularly  true  in  TMA  work  where  practical  limitations  preclude  "exact" 
modelling  of  the  actual  target  dynamics. 

Various  methods  have  been  suggested  to  prevent  the  growth  of  modelling 
errors;  fixed  memory  filters,  fading  memoiy  filters,  and/or  injecting  artificial 
plant  noise  into  the  system  are  but  a  few.  A  common  feature  of  these  methods 
is  that  they  all  attempt  to  prevent  divergence  by  discounting  the  influence  of 
past  data.  While  this  effectively  Increases  filter  response,  it  also  makes  the 
filter  more  susceptible  to  random  errors.  Consequently,  to  prevent  unnec¬ 
essary  sacrificing  of  noise  performance,  adaptive  control  is  required. 

The  need  for  adaptively  controlled  filters  has  provided  impetus  for  the 
development  of  numerous  "modified"  Kalman  filtering  schemes.  I-4  In  all  these 
schemes  the  adaptive  control  mechanism  is  a  function  of  the  Kalman  filter  data 
measurement  residuals.  This  might  be  expected  since  data  residuals  provide 
the  only  consistently  reliable  basis  for  assessing  solution  quality  and  detecting 
the  onset  of  divergence.  However,  it  is  unfortunate  that  the  importance  and 
practical  utility  of  residuals,  particularly  in  the  context  of  adaptive  filtering, 
have  not  been  adequately  stressed.  In  fact,  there  is  a  paucity  of  technical 
literature  devoted  to  this  subject. 

This  report  analyzes  the  Kalman  filter  and  its  associated  data  measure¬ 
ment  residuals.  By  utilizing  the  classical  method  of  least  squares,  a  deriva¬ 
tion  of  the  well  known  filter  equations  is  presented.  Although  this  particular 
approach  is  somewhat  unorthodox,  it  is  nevertheless  appealing  since  ancillary 
mathematical  formulas  that  are  later  needed  for  analyzing  the  residuals  emerge 
naturally.  Following  this,  an  important  relationship  between  the  system  per¬ 
formance  index  and  the  riata  measurement  residuals  is  developed  and  then  used 
to  demonstrate  the  suitability  of  residuals  for  providing  adaptive  control  and 
on-line  assessment  of  solution  quality.  In  addition,  pertinent  statistical  prop¬ 
erties  of  the  residuals  and  performance  index  are  deduced  and  utilized  to  es¬ 
tablish  a  basis  for  formulating  practical  adaptive  control  criteria.  Finally,  a 
simulation  example  which  demonstrates  filter  divergence  (e.g. ,  tracking  of  a 
maneuvering  target)  is  presented. 
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DERIVATION  OF  THE  KALMAN  FILTER  EQUATIONS 
DESCRIPTION  OF  THE  SYSTEM 

Consider  a  physical  system  which  can  be  described  mathematically  by  a 
set  of  linear  stochastic  vector  difference  equations  of  the  form 

X(k+1)  =  A(k+l,k)X(k)  +  B(k)W(k),  (1-a) 


Z(k+1)  =  H(k+l)X(k+l)  +  V(k+1),  k  =  0, 1.2,3 . 

where 
X(k) 


(1-b) 


-  p-dimensional  vector  which  statistically  describes  the  system 
states  at  time  Tk, 


A(k+l,k)  -  (pxp)  deterministic  one-step  transits 


B(k) 

W(k) 


-  (pxp)  deterministic  matrix  of  known  form, 

-  q-dimensional  vector  of  stochastic  inputs, 


tor  the  system. 


Z^k+1)  -  m-dimensional  vector  of  noisy  data  measurements  taken  at 

time  T, 


lk+l’ 


H(k+1) 

V(k+1) 


-  (mxp)  deterministic  measurement  matrix, 

-  m-dimensional  vector  of  additive  measurement  noise. 
The  pertinent  statistics  associated  with  this  system  are  as  follows: 

eJx(°)|  = 

Ejw(k)j  =  W(k,0) 
e|  V(k)  |  =  0 

eJ(X(0)  -  X(0,0)1[X(0)  -  X(0,0)’j  =  P(0) 

Ej(W(j)  -  W(j,0)l(W(k)  -  W(k,0)l'j  =6jkQ(k) 

Ejv(j)v'(k)j  =6jkR(k) 


(2-a) 
(2-b) 
(2-c) 
(2— d) 
(2-e) 
(2-f) 
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j  -  TE  4602 

•  v 


i. 


eJ[X<0)  -  X(0,0))lW(k)  -  W(k,0)l'{  =  0 

(2-g) 

e|(X(0)  -  X(0,0)]V*(k)  I  =  o 

(2— h) 

E  j  (W(J)  -  W(J,  0)]v'(k)  (  =  0 

(2-i) 

where  (l  j  =  k 

6jk=(o  j*k 

(2-J) 

and  e) •  (  denotes  the  statistical  expectation  operator. 


As  is  well  known  (e.  g. , 
given  by 


k-1 


Sage5),  the  formal  solution  to  equation  (1-a)  Is 


X(k)  -  A(k,  0)X(0)  +  £  A(k,  j+l)B(j)W(j),  k  =  1,2 .  (3) 


which  readily  reveals  the  statistical  nature  of  X(k)  through  its  functional 
dependence  on  the  random  vectors  X(0),W(0),W(1), . .  .W(k-l).  Consequently, 
estimates  of  X(k)  are  all  that  can  be  determined.  In  the  absence  of  additional 
information  it  is  evident  that  the  optimal  estimate  of  X(k)  is  Its  a  priori 
mean  value,  defined  by 


Ejx(k)  j  =  X(k,0)  =  A(k,0)X(0,0)  +  £  A(k,  j+l)B(J)W(j,0).  (4) 

1=0 


It  should  be  noted  that  this  equation  provides  an  a  priori  description  of  ex¬ 
pected  system  behavior.  Suppose,  however,  that  further  information  about 
the  system  states  is  provided  in  the  form  of  discrete  data  measurements 
Z(k)  which  are  linearly  related  to  X(k)  by  equation  (1-b).  If  such  information 
Is  properly  utilized,  it  is  reasonable  to  expect  that  better  estimates  of  X(k) 
will  result.  Numerous  techniques  exist  for  extracting  such  estimates  from 
measured  data.  One  of  the  most  powerful  i.<  the  Kalman  filter  which  will  now 
be  derived. 
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DEVELOPMENT  OF  THE  FILTER  EQUATIONS 

To  analyze  the  system  previously  described,  assume  that  n  data  measure¬ 
ments  Z(1),Z(2), . . .  Z(n)  are  available.  The  problem  at  hand  is  to  determine 
an  optimal  estimate  of  the  current  state  vector  X(n)  from  these  measure¬ 
ments  and  the  a  priori  information.  Here,  the  criteria  for  optimality  will  be 
defined  so  that  the  resulting  estimates  best  fit  the  measured  data  while  simul¬ 
taneously  minimizing  deviations  from  a  priori  expected  system  behavior. 

For  convenience,  the  following  notation  will  be  employed: 

X(k,n)  =  optimal  estimate  of  X(k)  based  on  n  data  measurements, 

W(k,n)  =  optimal  estimate  of  W(k)  based  on  n  data  measurements. 

It  can  be  seen  that  equation (4),  which  defines  X  (k,0),is  consistent  with  this 
notation  since  the  a  priori  mean  value  of  X(k)  represents  the  optimal  esti¬ 
mate  of  this  vector  when  no  data  measurements  are  available.  Finally,  to 
ensure  that  all  estimates  satisfy  the  general  system  equations,  the  following 
constraints  are  imposed: 

X(k+1 ,  n)  =  A(k+1 ,  k)X(k,  n)  +  B(k)W(k,  n)  ,  k  =  0, 1 , 2 , . . .  (n-1)  .  (5) 

To  complete  a  mathematical  formulation  of  the  problem,  it  is  necessary  to 
define  a  system  performance  index  J(n).  To  this  end,  let 


J(n)  =  l/2[X(0,n)  -  X(0, 0)l’p"1(0)[X(0,n)  -  X(0,0)J 
n-l  ,  . 

+  1/2  £  (W(l,n)  -  W(1,0)J  Q'NllfWd.n)  -  W(1,0)] 

1=0 

n  , 

+  1/2  £  [Z(l)  -  H(l)X(l,n)J  R  ‘(WZd)  -  H(l)X(l,n)J  .  (6) 

1=1 

It  might  be  noted  that  J(n)  is  a  weighted  quality  measure  of  the  estimation 
process  based  upon  the  optimality  criteria  employed.  The  first  two  terms  in 
equation  (6)  account  for  deviations  from  expected  system  behavior,  and  the 
last  term  accounts  for  cumulative  data-fit  errors.  Because  of  the  way  J(n) 
is  defined,  it  is  reasonable  to  expect  that  this  quantity,  or  some  function  thereof, 
will  provide  a  suitable  basis  for  assessing  solution  quality  on-line.  This  point 
is  explored  further  in  the  next  section. 
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An  optimal  estimate  X(n,  n)  of  the  current  state  vector  X(n)  may  now  be 
obtained  by  minimizing  J(n)  subject  to  the  constraints  imposed  by  equation  (5). 
However,  if  equation  (3)  is  utilized,  these  constraints  may  be  rewritten  in  the 
equivalent  format 

X(k,n)  =  A(k,0)X(0,n)  +  E  A(k,J+l)B(j)W(j,n),  k  =  l,2,...n,  (7) 

]=0 


which  reveals  that  J(n)  is  a  function  of  n+1  independent  vectors  X(0,n),W(0,n), 
W(l,n), . .  .W(n-l,n).  Consequently,  the  minimization  of  J(n)  will  require 
that 


5X(0,n) 


=  0, 


(8-a) 


0J(n) 

0W(k,n) 


=  0, 


k  =  0,1,2, . . . (n-1) . 


(8-b) 


Performing  the  operations  indicated  by  equations  (8)  yields 


P_1(0)[X(0,n)  -  X(0,0)]  =  E  A'(l,0)H'a)R'1a)IZ(l)  -  H(l)X(l,n)]  ,  (9-a) 

1=1 

n— i  * 

Q_1(k)[W(k,n)  -  W(k,0)J  =  {Z(\ )  -  Ha>xa,n)l, 

k  =  0,1,2,... (n-1)  .  (9-b) 

By  taking  the  transpose  of  equation  (7),  replacing  the  index  k  with  1,  and 
differentiating  the  result  with  respect  to  W(k,n),  it  follows  that 

OX(l.n)  i  o  l<(k+l) 

bW(k.n)  =  j  B'(k)A’(l,k+l)  U(k+1)  *  (10) 

Hence,  equation  (9-b)  reduces  to 

Q_1(k)[W(k,n)  -  W(k,0)]  =  B(k)  E  A'a.k+DHflJR'Vuza)  -  H(l)X(l,n)l  , 

l=k+l 

k  =  0, 1,2, . . .  (n-1) ,  (11) 
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The  determination  of  X(n,n)  first  requires  that  equations  (9-a)  and  (11), 
or  their  equivalent,  be  solved  for  the  unknown  vectors  X(0»n),  W(0,n),W(l,n), 
W(2,n),. .  ,W(n-l,n).  To  this  end,  it  proves  convenient  to  introduce  an  auxiliary 
vector  S(k,  n)  according  to  the  formula 

B(k)A(k, k+l)S(k, n)  =  Q_1(k)lW(k,n)  -  W(k,0)]  ,  k  =  0,1,2, ...  (n-1) .  (12) 

Substituting  this  expression  into  equation  (11)  and  utilizing  the  transition  ma¬ 
trix  property 

A(l,k+1)  =  A(k,k+l)A(l,k)  (13) 

leads  to  the  relationship 

S(k,n)  =  £  A  (1, k)H  a)R  (1)[Z(1)  -  H(l)X(l,n)],  k  =  0,1,2, . . .  (n-1).  (14) 
l=k+l 

If  the  index  k  in  equation  (14)  is  set  equal  to  zero  and  the  result  compared  to 
equation  (9-a),  the  following  boundary  condition  may  be  derived: 

S(0,n)  =  P_1(0)[X(0,n)  -  X(0, 0)).  (15) 

Similarly,  setting  k  =  (n-1)  in  equation  (14)  yields 

S(n-l,n)  =  A  (n,n-l)H  (n)R_1(n)[Z(n)  -  H(n)X(n,n)] .  (16) 

Again  referring  to  equation  (14),  it  can  be  seen  that 

S(k+1 , n)  =D  A'd.k+lJHajR'^UtZd)  -  H(l)X(l,n)l 

l=k+2 

=  A'(k,k+l)S(k,n)  -  H’(k+l)R"1(k+l)(Z(k+l)  -  H(k+l)X(k+l,n)l, 

k  =  0,1,2,...  (n-1),  (17) 

where  the  transition  matrix  identity 

A(k,k)  =  I ,  k  =  0,1,2,...,  (18) 
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has  been  utilized.  Note  that  for  k  =  (n— 1)  a  new  vector  S(n,n)  now  appears 
in  equation  (17)  and  is  defined  by 

S(n,n)  =  A(n-l,n)S(n-l,n)  -  H'(n)R_1(n)lZ(n)  -  H(n)X(n,n)].  (19) 


However,  substituting  equation  (16)  into  (19)  and  noting  that 

A(l,  k)  =  A-1(k,l)  (20) 

leads  to  the  boundary  condition 

S(n,n)  =  0  .  (21) 

Equations  (7),  (12),  (IS),  (17),  and  (21)  comprise  a  coupled  set  of  linear  inhomo¬ 
geneous  difference  equations  and  associated  boundary  conditions  that  are  math¬ 
ematically  equivalent  to  equations  (9-a)  and  (11).  .This  set,  which  is  summarized 
below,  may  therefore  be  used  to  determine  X(n,n). 

X(k+1 , n)  =  A(k+l,k)X(k,n)  +  B(k)W(k,0)  +  B(k)Q(k)B'(k)A’(k.k+l)S(k,n) ,  (22-a) 

S(k+1 , n)  =  A'(k,k+l)S(k,n)  -  H'(k+l)R"1(k+l)[Z(k+l)  -  H(k+l)X(k+l,n)] ,  (22-b) 

S(0,n)  =  P_1(0)[X(0,n)  -  X<0,0)]  ,  (22-c) 

S(n,n)  =0,  k  =  0,1,2, ...(n-1).  (22-d) 


It  is  interesting  to  note  that  the  above  equations  are  similar  to  the  discrete 
canonical  equations  derived  by  Sage5  using  dynamic  optimization  theory. 

Both  sets  describe  linear  two-point  boundary  value  problems  which  may  be  re¬ 
solved  by  the  method  of  discrete  invariant  imbedding.  5  However,  an  alter¬ 
nate  procedure  will  be  employed  here,  which  yields  the  desired  vector  X(n>n) 
in  a  relatively  simple  fashion.  Thus,  assume  that  one  additional  data  measure¬ 
ment  is  taken  at  time  Tn+i .  The  data  measurement  sequence  is  now  defined  by 
{ Z(l),  55(2), . . .  Z(n+1) }  and  equation  set  (22)  is  replaced  by 


X(k+l,n+l)  =  A(k+l,k)X(k,n+l)  +  B(k)W(k,0)  +  B(k)Q(k)B’(k)A(k,k+l)S(k,n+l), 

(23-a) 

S(k+l,n+l)  =  A  (k, k+l)S(k,n+l)  -  H'(k+l)R_1(k+l)[Z(k+l)  -  H(k+l)X(k+l,n+l)], 


S(0, n+1)  =  P"1(0)fX(0,n+l)  -  X(0,0)J, 


(23-b) 

(23-c) 


k  =  0, 1,2, . .  ,n. 


S(n+l,n+l)  =  0  , 


(23-d) 
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Next,  define 


AX(k)  =  X(k,  n+1)  -  X(k,  n) ,  (24-a) 

AS(k)  =  S(k,n+1)  -  S(k,n) .  (24-b) 

Subtracting  the  expressions  in  equation  set  (23)  from  the  corresponding  expres¬ 
sions  in  equation  set  (22)  and  utilizing  equations  (24)  then  yields 

AX(k+l)  =  A(kH,k)AX(k)  +  B(k)Q(k)B'(k)A(k.k+l)AS(k) ,  (25-a) 

AS(k+l)  =  A(k,k+l)AS(k)  +  H'(k+l)R"l(k+l)H(k+l)AX(k+l) ,  (25-b) 

AS(0)  =  P-1(0)AX(0) ,  (25-c) 

AS(n)  =  S(n,n+l).  (25-d) 

Though  not  immediately  evident,  the  homogeneous  equation  set  (25)  is 
much  simpler  to  solve  than  the  original  equation  set  (22).  To  illustrate,  assume 
a  solution  of  the  form 

AX(k)  =  P(k)AS(k) ,  k  =  0,1,2,...  (n),  (26) 


where  P(0)  is  known,  and  P(l),  P(2), . . .  P(n)  are  matrices  to  be  determined. 
(Note  that  a  similar  type  homogeneous  solution  of  the  form  X(k,n)  -  P(k)S(k,n) 
cannot  be  applied  to  equation  set  (22).)  Subitituting  equation  (26)  in  (25-a)  and 
(25-b), respectively,  and  utilizing  equation  (20)  then  yields 

P(k+l)AS(k+l)  =  N(k+l)A(k,k+l)AS(k) ,  (27-a) 

(P-1(k+l)  -  H'(k+l)R*1(k+l)H(k+l)lP(k+l)AS(k+l)  =  A'(k,k+l)AS(k),  (27-b) 

where 

N(k+1)  =  A (k+1 ,  k)  P(k)  A  (k+l ,  k)  +  B(k)Q(k)B'(k).  (28) 

In  order  to  satisfy  equation  set  (27)  for  all  allowable  values  of  k ,  it  is 
necessary  that 


P(k+1)  =  (N(k+1)  +  H,(k+l)R"1(k+l)H(k+l)r1.  (29) 

In  addition,  equations  (24-a),  (25-d),  and  (26)  may  be  combined  to  produce  the 
expression 


X(n,n+1)  =  X(n,n)  +  P(n)S(n,  n+1) . 


(30) 


8 


TR  4502 


Equation  (23-a),  coupled  with  equations  (28)  and  (30),  then  yields  the  relation¬ 
ship 

X(n+l,n+l)  =  A(n+l,n)X(n,n)  +  B(n)W(n,0)  +  N(n+1)A  (n,n+l)S(n,n+l).  (31) 

Note  also  that  equations  (23-b)  and  (23-d)  lead  to 

A'(n,n+l)S(n,n+l)  =  H  (n+l)R“1(n+l)[Z(n+l)  -  H(n+l)X(n+l,n+l)].  (32) 

Finally,  substituting  equation  (32)  into  equation  (31)  to  eliminate  S(n,n+1)  and 
performing  some  algebraic  manipulation  produce  the  desired  results: 

X(n+l,n+l)  =  X(n+l,n)  +  P(n+l)H'(n+l)R”1(n+l)IZ(n+l) 

-  H(n+l)X(n+l,n)l,  (33-a) 


Xfn+l,n)  =  A (n+1 , n)X(n, n)  +  B(n)W(n, 0)  .  (33-b) 

Equations  (28),  (29),  and  (33)  provide  a  mathematical  procedure  for  opti¬ 
mally  estimating  the  current  state  of  a  system  from  measured  data  and  a  priori 
statistical  information.  As  expected,  these  equations  are  equivalent  to  the 
Kalman  filtering  formulas  and  can  be  put  into  the  standard  format  by  utilizing 
the  following  matrix  identities: 

G(n+1)  =  N  (n+l)H  *  (n+1)  [H(n+1  )N  (n+1  )H '  (n+1)  +  R(n+l)f 1 

*  P(n+1)H  (n+l)R_1(n+l),  (34-a) 

P(n+1)  =  [N-1(n+l)  +  H'(n+l)R"1(n+l)H(n+l)l"1 

=  [I  -  G(n+l)H(n+l)]N(n+l),  (34-b) 


For  completeness,  the  Kalman  filter  equations  are  summarized  on  the  following 
page  in  algorithm  form. 
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1  Kalman  Filter  Equations 

_ ^,0)  -  initial  estimate  of  state-vector 

P(&  -  initial  estimate  of  state-vector  covariance  matrix 


X(n+1, n)  =  A(n+l,n)X(n, n)  +  B(n)W(n,0) ,  (35-a) 

N(n+1)  =  A(n+l,n)P(n)A  (n+l,n)  +  B(n)Q(n)B  (n)  ,  (35-b) 

G(n+1)  =  |tf(n+l)H\n+l)[H(n+l)N>H)H'<n+l)  +  R(n+l)]_1,  (35-c) 

X(n+l,n+L)  =  X(n+l,n)  +  G(n+l)(Z(n+l)  -  H(n+l)X(n+l,n)]  ,  (35-d) 

P(n+1)  =  [I  -  G(n+l)H(n+l))N(n+l) ,  n  =  0,1.2 .  (35-e) 
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ANALYSIS  OF  THE  KALMAN  FILTER 
DATA  MEASUREMENT  RESIDUALS 


REPRESENTATION  OF  THE  PERFORMANCE  INDEX 
IN  TERMS  OF  RESIDUALS 

Data  measurement  residuals,  which  arise  naturally  during  computations 
performed  by  the  Kalman  filter,  are  defined  by  the  equation 


Y(k)  =  Z(k)  -  H(k)X(k,  k-1) ,  k  -  1, 2, ... .  (36) 

In  effect,  Y(k)  is  a  measure  of  the  error  between  the  actual  data  measure¬ 
ment  vector  at  time  Tjj  and  the  best  available  prediction  of  that  vector.  When 
taken  collectively, these  residuals  provide  reliable  indication  of  how  well  the 
state-vector  estimates  "fit"  the  measured  data.  In  turn,  such  information  can 
be  effectively  utilized  to  assess  solution  quality.  However,  it  is  important  to 
recognize  that  since  the  residuals  are  random  quantities  they  only  provide  in¬ 
formation  of  a  statistical  nature.  The  magnitude  of  any  one  residual  has  little 
significance  by  itself. 

In  the  preceding  section,  it  was  pointed  out  that  the  performance  index  J(n) 
provides  a  basis  for  assessing  solution  quality  on-line.  Since  the  residuals 
may  also  be  used  for  this  purpose,  it  would  appear  that  these  quantities  are  re¬ 
lated  to  J(n)  in  some  way.  To  determine  this  relationship,  first  note  that  equa¬ 
tions  (€',  and  (9-a)  may  be  combined  to  yield  the  expression 


n-1 


J(n)  =  1/2  £  QWM)  -  W0.°>] VV)[WU,n)  _  W(1,0)] 

1=0 


n 

+  \nE  IZ<1)  -  H(1)X(1,  n)  «■  H(1)A(1, 0)X(0, n) 
1=1  “ 


- H(i)A(l, 0)X(0,  0)]'r_1(1)[Z(1)  - H(l)X(l,n)J. 


However,  equations  (4)  and  (7)  reveal  that 


(37) 


X(k,n)  -  A(k,  0)X(0,n)  +  A(k,0)X(0,0) 


k-1 

=  X(k,0)+L  A(k,j+l)B(j)[W(J,n)  -  W(J,0)].  (38) 

j=0 
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Consequently,  J(n)  may  be  rewritten  in  the  form 

J(n)  -  1/2  £1|Wa.“)-W(1.0)l,Q’1a)(W(l,n)-W(l,0)] 
1-0 


n  1-1 


i  i  i 


+1/2 E  E  iwa.nj-wa.oi'B'oiA’a.j+^HaiR^aitzaj-Hdixa.n)] 

1=11-0 


+i/2 £  [za>  -  Haijca.o)]  r^auza)  -  H<i)X(i,n)i .  m 

i-i 

The  eeoond  term  on  the  right-hand  aide  of  equation  (39)  can  be  simplified  by 
noting  that 

£  E  tWa,n)-Wa.O)l,B,(j)A,(l,J+l)H,a)R’1fl)lZ(l)-Hfl)Xa.n)l 

1=1  J-0 

-  E1  Kfl.n)  -  Wfl.  on'Bd)  1}  a'o.I+Dh'oiR^OIIZU)  -  H(])X(]»n>] 

1-0  J-l+1 

-  W'*)  -  Wfl.  0)l  Q-l(l)(Wfl#n)  -  Wfl,  0)1  ,  (40) 

1-0 

where  the  final  result  derives  from  equation  (11).  If  equation  (40)  is  now  sub¬ 
stituted  Into  equation  (39)  It  follows  that 

n  »  i 

J(n)  -  1/2  E  HO)  -  « (1)  Xfl.  0)1  R'A(1)(Z(1)  -  H(l)X(l,n)l .  (41) 

1-1 

Next,  note  that  equations  (35-a),  (35-d),  and  (36)  may  be  oomblned  to 
produce  the  expression 

X(k+l,k+l)  -  A(k+l,k)X(k,k)  +  B(k)W(k, 0)  +  G(k+l)Y(k+l) ,  (42) 


whioh  can  also  be  written  In  the  equivalent  form 

k-1  k 

X(k,k)  -  A(k,  0)X(0, 0)  +  E  A(k,  J+1)B(J)W(J,  0)  +  E  A(k,J)G(J)Y(J) .  (43) 

J-0  J-l 
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Comparing  this  equation  with  equation  (4)  reveals  that 

k 

X(k,  0)  =  X(k,  k)  -£  A(k,  J)G(j)Y(J) .  (44) 

1=0 


Finally,  substituting  equation  (35-d)  into  equation  (44)  to  eliminate  X(k»k) 
yields 


k 

X(k,  0)  =  X{k,k-1)  +  G(k)Y(k)  -  £  A(k,  j)G(J)Y(j)  ,  (45) 

1=1 


from  which  it  is  readily  deduced  that 

k 

Z(k)  -  H(k)X(k,  0)  =  [I  -  H(k)G(k)]Y(k)  +  H(k)  £  A(k,  j)G(j)Y(j)  .  (46) 

1=1 


If  equation  (46)  is  now  substituted  into  equation  (41),  the  performance  Index 
J(n)  takes  the  form 

'  »  •  _i 

J(n)  *  1/2  £  Y  (1)P  -  G  (1)H  (1)]R  (1)[Z(1)  -  H(l)X(l,n)] 

1-1 


f  t  f 


+  1/2  IL  Y  ( j  )G  (J)  A  (1 ,  j)  H '  (1)R  ~ 1  (1)  [Z^(l)  -  H(l)X(l,n)]  .  (47) 

1=1  j=l 


This  cumbersome  expression  can  be  greatly  simplified  by  noting  that 


£  £  Y'aG'ojA'a.DH'aiR^aiizd)  -  H(i)xa,n)] 

1=1  j=l 

=  £  y'(1)g’(1)L  A'o.DH’dJR'^JllZd)  -H(j)X(l,n)l 

1=1  j=l 


=  £  Y^DG'dijH’dJR'VllZd)  -  H(l)Xd.n)]  +  S(l,n)  j  ,  (48) 

1=1 

where  the  final  results  derive  from  equations  (14)  and  (18).  Hence,  by  com¬ 
bining  equations  (47)  and (43),  the  expression  for  J(n)  reduces  to 
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l 


J(n)  =  1/2  E  l'a) I R'1  a) (Z(l)  -  H(l)X(l,n)]  +  G  (l)S(l,n)|  • 
1=1 


(49) 


However,  since  n  is  an  arbitrary  Integer,  it  follows  Immediately  that 


Subtracting  equation  (49)  from  equation  (50)  then  yields  the  recursive  formula 

2-,  »  i  _i 

J(n+1)  -  J(n)  =  1/2  El  (1)(C  (l)AS(l)  -  R  (1)H(1)AX(1)] 

1=1 

+  l/2Y'(n+l)jR_1(n+l)[Z(n+l)  -  H(n+l)X(n+l,n+l)]  +  G(n+l)S(n+l,n+l) j  ,  (51) 

where  AX(1)  and  AS(1)  are  defined  by  equation  set  (24).  The  next  step  in  the 
reduction  process  involves  substituting  equations  (25-d)  and  (26)  into  equation 
(51)  to  obtain 

"  . 

J(n+1)  -  J(n)  =  1/2  E  I  (1)(G  (1)  -  R  (1)H(1)P(1)]AS<1) 

1=1 

+  1  /2  Y* (n+1  )R_1  (n+1 ) (Z (n+1)  -  H(n+l)X(n+l,n+l))  .  (52) 

Note,  however,  that  P(l)  and  R(l)  are  both  symmetric  matrices;  consequently, 
the  identity 


g'oc)  -  R_1(k)H(k)P(k)  =  0  ,  k=  0,1,2,...,  (53) 

is  easily  derived  by  taking  the  matrix  transpose  of  equation  (34-a).  This 
relationship  shows  that  the  series  appearing  in  equation  (51)  vanishes.  In  addi¬ 
tion,  the  last  remaining  term  on  the  right-hand  side  of  equation  (51)  can  be  sim¬ 
plified  by  utilizing  equations  (35-c),  (35-d),  (35-e),  and  (36).  The  result  is 

Y  (n+l)R-1(n+l)(Z(n+l)  -  H(n+l)X(n+l,n+l)] 

=  Y* (n+1) R'1  (n+1) [I  -H(n+l)G(n+l)]Y(n+l) 

=  Y'(n+l)(H(n+l)N(n+l)H  (n+l)  +  R(n+l))_1Y(n+l)  .  (54) 
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The  preceding  equations  may  now  be  combined  to  produce  a  simple  recursive 
formula  for  J(n);  that  is 

J(n+1)  =  J(n)  +  l/2Y(n+l)[H(n+l)N(n+l)H(n+l)  +  R(n+l)f  JY(n+l)  , 

n  =  0,1,2, ... .  (55) 

An  initial  condition  for  J(n)  can  also  be  deduced  from  equation  (49).  More 
precisely,  since 

J(l)  =  l^Y’dljR'^DlZd)  -  H(1)X(1, 1)]  +  G(1)S(1,1)| 

=  1/2Y  (1)R-1(1)(Z^(1)  -  H(1)X(1, 1)1 
=  1/2Y,(1)R“1(1)II  -  H(1)G(1))Y(1) 

=  l^YdimdjNdjH’d)  +  R<i)f  ^u) , 
it  follows  from  equation  (55)  that 

J(0)  =  0  .  (56) 

Although  equations  (55)  and  (56)  are  convenient  for  numerical  work,  a  closed 
form  representation  of  J(n)  is  often  more  desirable  for  analytical  purposes. 

Such  a  representation  takes  the  form 

i  *  _i 

J(n)  =  1/2  £  Y  (1)[H(1)N(1)H  (1)  +  R(l)]  AY(1)  .  (57) 

1=1 


A  comparison  of  this  expression  with  equation  (6)  shows  that  the  data 
measurement  residuals  also  provide  a  weighted  quality  measure  of  the  estima¬ 
tion  process  based  upon  the  optimality  criteria  employed.  As  such,  it  now 
becomes  evident  why  these  residuals  can  be  utilized  effectively  to  assess  filter 
performance. 

An  important  feature  of  equation  (57)  not  to  be  overlooked  is  that  all  quan¬ 
tities  needed  to  calculate  J(n)  can  be  extracted  directly  from  the  Kalman  filter 
equations  without  further  computation.  In  addition,  once  the  terms  in  the  series 
have  been  evaluated  they  need  not  be  recomputed  with  each  new  data  measure¬ 
ment.  Consequently,  for  practical  applications,  the  utilization  of  equation  (57)  — 
in  lieu  of  equation  (6)  —  can  result  in  a  significant  saving  of  computer  time. 
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STATISTICAL  PROPERTIES  OF  THE  RESIDUALS 

Earlier  In  the  discussion,  it  was  noted  that  the  data  measurement  residuals 
are  random  variables  and  can  only  provide  information  of  a  statistical  nature. 
The  performance  index  also  possesses  these  characteristics.  Consequently, 
knowledge  of  the  statistical  properties  of  both  J(n)  and  Y(k)  is  required 
before  meaningful  information  can  be  extracted  from  these  quantities. 

The  mean  value  of  Y(k)  can  be  determined  rather  easily  by  combining 
equations  (1-b)  and  (46)  to  obtain 


k 

H(k)(X(k)  -  X(M>1  +  V(k)  =  [I  -  H(k)G(k)]Y(k)  +  H(k)  £  A(k,J)G(J)Y(J)  .  (58) 

J=1 


Taking  the  expected  value  of  this  equation  and  utilizing  equations  (2-c)  and  (4) 
then  yields 


fl  -  H(k)G(k)lEjY(k)j  +  H(k)  £  A(k,J)G(j)E|Y(j)  j  =  0  . 


(59) 


Since  the  index  k  in  equation  (59)  Is  arbitrary,  it  can  easily  be  shown  by 
induction  that 


E|Y(k)f=0,  k  =  1,2 .  (60) 

In  order  to  determine  the  covariance  of  Y(k),  it  is  first  necessary  to  prove 
the  following  identity: 

E |  [X(k)  -  X(k,k-l)](X(k)  -  X<k,k-l)l'|  =  N(k) .  (61) 

To  this  end,  assume  that  equation  (61)  is  true  for  some  fixed  value  of  k.  By 
combining  equations  (1-a),  (1-b),  (35-a),  and  (35-d)  it  can  also  be  shown  that 

X(k+1)  -  X(k+l,k)  =  A(k+l,k)[X(k)  -  X(k,k>]  +  B(k)[W(k)  -  W(k,0)) 

=  A(k+l,k)(X(k)  -  X(k,k-1)  -  G(k)Z(k)  +  G(k)H(k)X(k,k-l)] 

+  B(k)(W(k)  -  W(k,0)J 
=  A(k+l,k)[I  -  G(k)H(k)](X(k)  -  X(k.k-l)] 

-  A (k+1 ,  k)G (k)V(k)  +  B(k)[W(k)  -  W(k,  0)] .  (62) 
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Consequently,  the  statistical  properties  depicted  in  equation  set  (2)  and  equa¬ 
tion  (61)  can  be  utilized,  together  with  equation  (62),  to  produce  the  expression 

E  j  (X(k+1 )  -  X(k+1 ,  k)  ](X(k+l)  -  X(k+1 ,  k) ) '  | 

=  A(k+l,k)|[I  -  G-(k)H(k)]N(k)[I  -  G(k)H(k)]'  +  G(k)R(k)G'(k)j  A'(k+l,k) 

+  B(k)Q(k)B*(k)  .  (63) 


Note,  however,  that 

P(k)  =  [I  -  G(k)H(k)]N(k) 

=  [I  -  G(k)H(k))N(k)[I  -  G(k)H(k)]'  +  G(k)R(k)G'(k)  .  (64) 

Consequently,  equation  (63)  reduces  to 

Ej[X(k+l)  -  X(k+l,k)UX(k+l)  -X(k+l,k))'| 

=  A  (k+1 ,  k)  P(k)  A ' (k+1 , k)  +  B(k)Q(k)B'(k)  =  N(k+1) ,  (65) 

where  the  final  result  derives  from  equation  (35-b).  It  has  therefore  been 
shown  that  if  equation  (61)  is  valid  for  any  value  of  k,  then  so  is  equation  (65). 
The  last  remaining  step  in  this  induction  proof  is  to  show  that  equation  (61)  is 
valid  for  some  value  of  k.  In  particular,  let  k  =  1.  For  this  case,  equations 
(1-a)  and  (35-a)  can  be  combined  to  yield 

X(l)  -  X(l, 0)  =  A(l, 0)(X(0)  X<0, 0)1  +  B(0)[£(0)  -  W(0, 0)]  .  (66) 


If  equation  set  (2)  is  again  utilized,  together  with  equations  (35-b)  and  (54),  it 
then  follows  that 

Ej[X<l)  -  X(l, 0)12(1)  -  SI, 0)1 ' | 

=  A(1,0)P(0)a'(1,0)  +  B(0)Q(0)^'(0)  =  N(l) ,  (67) 

and  the  validity  of  equation  (61)  is  established. 
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The  covariance  of  Y(k)  can  now  be  determined  in  a  simple  manner  by 
recalling  that 

Y(k)  =  H(k)[X(k)  -  X(k,  k-1)]  +  V  (k)  .  (68) 

As  such,  the  desired  result 

Ej  Y(k)Y'(k)|  =  H(k)N(k)H(k)  +  R(k)  (69) 

can  be  deduced  Immediately  from  equation  set  (2)  and  equation  (61).  It 
might  also  be  noted  that 

EjX(k)(H(k)N(k)H(k)  +  R(k)]_1Y(k)  j 

«  Trace  J  lH(k)N(k)H'(k)  +  R(k)]_1E{  Y(k)Y'(k)|j 

=  Trace  j  [H(k)N(k)H'(k)  +  R(k))"1[H(k)N(k)H'(k)  +  R(k)]  [ 

-Trace  |lf  =  m  ,  (70) 

where  m  is  the  dimension  of  the  measurement  vector  Z(k).  With  the  aid  of 
equation  (70),  the  mean  value  of  the  performs  ace  index  can  be  easily  determined 
as  follows: 

^|X(DW1)N(1)H(1)  +  R(l)l""1Y(l)  ( 
m  =  (1/2)  mn  .  (71) 


Mn 

=  l/2l 
1=1 

=  1/2^ 
1=1 


A  final  statistic  of  importance  is  the  variance  of  J(n).  In  general,  a 
closed  form  expression  for  this  quantity  cannot  readily  be  obtained  since  fourth 
order  statistical  moments  must  be  computed.  However,  if  all  statistics  are 
Gaussian,  then  J(n)  will  be  Chi-Square  distributed6  with  m«n  degrees  of 
freedom.  Under  these  conditions,  the  variance  of  J(n)  is  given  by 

E  j  [J(n)  -  E  { J(n)  >  ]2  j  =  (l/2)mn  .  (72) 
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It  Is  important  to  note  that  this  statistic  provides  a  measure  of  the  permissible 
deviation  in  J(n)  from  its  expected  value.  Unfortunately,  since  larger  and 
larger  deviations  are  permitted  as  the  number  of  data  measurements  increases, 
it  becomes  progressively  more  difficult  to  extract  reliable  information  from 
J(n).  More  precisely,  after  many  measurements  have  been  processed,  large 
deviations  in  J(n)  cannot  be  attributed  solely  to  filter  divergence,  but  instead 
may  result  simply  from  the  accumulation  of  random  errors.  The  con¬ 
comitant  ambiguity  makes  it  impossible  to  accurately  assess  the  true  filter 
status. 

Fortunately,  there  are  various  ways  to  circumvent  this  difficulty.  The 
simplest  technique  is  to  introduce  a  "modified"  performance  index  of  the  form 

L(n)  =  1/2  Z  A1) Y'fl)lHa)N(l)H'(i)  +  R(1)]_1Y(1)  -  m|  ,  (73) 

1=1 


where  y  is  a  weighting  factor  which  satisfies  the  inequality 

0<y <  1 . 

The  mean  and  variance  of  L(n)  are  given  by 

E )  L(n)  (  =0  , 


(74) 

(75-a) 

(75-b) 


Note  that  these  statistics  remain  bounded  as  the  number  of  data  measurements 
increase.  In  fact 


lim  E  jL2(n)  (  = 


m 


n-»® 


2  * 
2(1-/  ) 


(76) 


For  finite  values  of  n,  any  desired  variance  in  the  range 


f<E)L2(n,j<m 


(77) 


can  be  obtained  by  choosing  y  appropriately.  This  flexibility  is  achieved  by 
exponentially  discounting  past  filter  performance  so  that  L(n)  is  influenced 
much  more  by  current  events  than  is  J(n).  As  such,  L(n)  will  provide  a  re- 
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liable  quantitative  measure  of  current  solution  quality  which  can  be  effectively 
utilized  for  adaptive  control.  For  example,  if  the  values  of  L(n)  are  statis¬ 
tically  consistent,  then  it  can  be  safely  assumed  that  the  Kalman  filter  is  per¬ 
forming  satisfactorily.  However,  if  L(n>  becomes  statistically  inconsistent 
during  the  course  of  operation  (e.g. ,  L(n)  continually  exceeds  its  one  sigma 
limit),  this  would  Indicate  that  the  filter  is  diverging,  and  appropriate  preventa¬ 
tive  action  should  be  taken.  It  might  be  noted  that  the  effectiveness  of  any 
divergence  prevention  scheme  can  also  be  easily  assessed  by  continually  moni¬ 
toring  L(n)  for  statistical  consistency. 

For  convenience,  a  recursive  algorithm  for  numerically  computing  L(n) 
as  well  as  a  summary  of  its  statistical  properties  is  presented  below: 


L(0)  =  0  , 


(78-a) 


Y(n+1)  =  Z(n+1)  -  H(n+l)X(n+l,n)  ,  (78-b) 

L(n+1)  =  YL(n)  +  l^Jy^n+lJIHin+lJN^+llH^n+l)  +  R(n+l)f  ^(n+l)-  m|,(78-c) 
EjY(n+l)j  =  0,  (78-d) 

Ej  Y(n+l)Y'(n+l)  j  =  H(n+l)N(n+l)H'(n+l)  +  R(n+1)  ,  (78-e) 

EjL(n+l)j  =  0  , 


(78— D 
(78-g) 


e|lVi)(  =  ®  0<y<1  * 

where  m  is  the  dimension  of  the  data  measurement  vector,  and  n  =  0,1,2,.. . 
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SIMULATION  RESULTS 

To  demonstrate  the  advantages  of  using  a  modified  performance  index  for 
adaptive  control,  a  target  motion  analysis  (TMA)  experiment  was  conducted. 
Since  target  maneuvers  are  a  primary  cause  of  divergence  in  present  TMA  sys¬ 
tems,  the  geometry  depicted  in  figure  1  was  chosen.  In  this  geometry,  target 
maneuvers  consisting  of  180  degree  turns  occur  at  approximately  6  minutes  and 
12  minutes.  Ownship  motion  is  as  shown  and  the  turning  rate  of  both  vehicles 
is  3°  per  second.  *  Ownship  speed  is  18  knots  and  target  speed  is  20  knots. 

Two  experiments  were  performed  with  this  geometry. 

First,  an  extended  non-adaptive  Kalman  filter-TMA  algorithm  which 
utilizes  a  constant  target  acceleration  plant  description  was  employed.  Data 
measurements  consisted  of  bearing  and  range  rate  information.  Figures  2,  3, 
and  4  are  plots  of  the  target  parameter  estimate  errors  versus  time.  In  figures 
5  and  6,  normalized  valuest  of  the  performance  index  J(n),  defined  by  equation 
(57),  and  the  modified  performance  index  L(n),  defined  by  equation  (73),  are 
presented.  From  examination  of  figures  2,  3,  and  4,  it  is  evident  that  the  esti¬ 
mation  errors  become  unacceptably  large  after  the  first  target  maneuver.  In 
addition,  deviations  (from  expected  behavior)  apparent  in  figures  5  and  6  indi¬ 
cate  that  the  filter  statistics  are  no  longer  consistent  with  their  a  priori  expected 
values.  Clearly,  the  target  maneuvers  have  caused  the  Kalman  filter  to  diverge. 

In  the  second  experiment,  an  extended  adaptive  Kalman  filter-TMA  algo¬ 
rithm  which  also  utilizes  a  constant  target  acceleration  plant  description  was 
employed.  Here,  however,  plant  noise  was  injected  into  the  system  to  accom¬ 
modate  mis-modelling  errors  introduced  by  target  maneuvers.  The  plant  noise 
covariance  matrix  is  adaptively  varied  according  to  the  behavior  of  the  modified 
performance  index.  Figures  7,  8,  and  9  are  plots  of  the  target  parameter 
estimate  errors  versus  time  for  this  case.  Figures  10  and  11  are  the  respec¬ 
tive  analogues  to  figures  5  and  6.  Inspection  of  figures  7,  8,  and  9  reveals 
that  the  error  divergence  observed  in  figures  2,  3,  and  4  is  no  longer  present; 
and  examination  of  figures  10  and  11  indicates  that  the  adaptive  control  mech¬ 
anism  has  effectively  realigned  the  filter  statistics  with  their  a  priori  expected 
values.  Further  examination  of  figures  10  and  11  also  reveals  the  sluggish 
nature  of  the  performance  index  J(n)  and  the  more  responsive  behavior  of  the 
modified  performance  index  L(n). 


♦Turning  radius  of  ownship  is  not  shown  in  figure  1. 

fin  this  simulation,  the  performance  indices  J(n)  and  L(n)  have  been 
normalized  so  that  both  have  expected  values  of  unity  for  all  n. 


21 


TR  4502 


Figure  1.  Target-Ownship  Geometry 
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Figure  7.  Target  Range  Error  vs  Time  (Adaptive  Kalman  Filter) 


Figure  8.  Target  Course  Error  vs  Time  (Adaptive  Kalman  Filter) 


Figure  9.  Target  Speed  Error  vs  Time  (Adaptive  Kalman  Filter) 
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Figure  10.  Performance  Index  vs  Time  (Adaptive  Kalman  Filter) 


Figure  11.  Modified  Performance  Index  vs  Time  (Adaptive  Kalman  Filter) 
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CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 

An  important  relationship  between  the  Kalman  filter  data  measurement 
residuals  and  the  associated  system  performance  index  has  been  established. 
The  computational  simplicity  of  the  relationship  not  only  permits  on-line  eval¬ 
uation  of  the  performance  index,  but  is  further  exploited  to  deduce  pertinent 
statistical  properties  of  the  index.  These  statistics  are  subsequently  utilized 
to  formulate  practical  adaptive  control  criteria  and  to  modify  the  performance 
index  to  enhance  its  reliability  as  a  solution  quality  indicator.  The  utility  of 
the  modified  performance  index  for  assessing  filter  status  and  for  adaptively 
regulating  the  plant  noise  covariance  matrix  has  been  discussed  and  demon¬ 
strated  via  laboratory  simulation. 

Even  though  the  simulation  results  are  preliminary  they  clearly  illustrate 
that  a  significant  improvement  in  filter  performance  can  be  achieved  when  this 
type  of  adaptive  control  mechanism  is  appended. 

Although  the  results  presented  here  are  encouraging,  much  work  remains 
to  be  done.  In  particular,  a  quantitative  functional  relationship  between  the 
modified  performance  index  and  the  plant  covariance  matrix  needs  to  be  estab¬ 
lished.  A  study  t°  determine  how  at-sea  data  affect  the  behavior  of  the  modi¬ 
fied  performance  index  should  also  be  conducted.  Finally,  techniques  to  adap¬ 
tively  compute  both  the  mean  and  covariance  of  the  plant  noise  should  be  further 
explored.  Current  efforts  are  being  directed  toward  the  solution  of  these 
problems  and  the  results  of  these  investigations  will  be  documented  in  future 
reports. 
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